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ABSTRACT 

A global Markov Chain Monte Carlo analysis of published eclipse photometry and 
radial velocities is presented for the transiting planet HD 80606b. Despite the lack of 
a complete transit light curve, the size of the planet is measured with a good level of 
precision {Rp = 1.04^g Qg Rjup), while the orbital parameters are refined. This global 
analysis reveals that the orbital axis of the planet is significatively inclined relative to 
the spin axis of the host star {(3 = — 59;1^28 "i^g), providing a compelling evidence that 
HD 80606b owes its peculiar orbit to the Kozai migration mechanism. 
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1 INTRODUCTION 

Among the known transiting gazeous planets, HD 80606b is 
definitely an unique case. Detected by radial velocity mea- 
surements (Naef et al. 2001), this massive planet (M ~ 4 
Mj) orbits in a very eccentric orbit (e ~ 0.93) around one 
of the two solar-type components of a wide binary (~ 1000 
AU). With P ~ 111 days, it is by far the known transiting 
planet with the largest period. Laughlin et al. (2009; herafter 
L09) observed this planet with Spitzer during its periastron 
passage and detected its occultation, implying an inclination 
close to 90°. A partial photometric transit was firmly de- 
tected a few months later by several ground-based telescopes 
(Fossey et al. 2009, hereafter F09; Garcia-Melendo & Mc- 
CuUough 2009, hereafter G09; Moutou et al. 2009, hereafter 
M09). This transit was also detected by M09 in spectroscopy 
through the observation of the Rossiter-McLaughlin effect 
(Queloz et al. 2000). 

AU three transit detection papers report a size close to 
Jupiter's for the planet. Still, the shape and duration of the 
transit are poorly constrained by the ground-based partial 
transit data. G09 considered a grazing transit as unlikely, 
basing on the similar transit depths measured in B-band 
(G09) and J?-band (F09, M09), but did not consider the 
possibility of an imperfect normalization or the presence of 
a trend or low-frequency noise (due, e.g., to an imperfect 
colored extinction correction or to a low-frequency variabil- 
ity of the star) in one or more of these ground-based time- 
series. Assessing the infiuence of such an effect is especially 
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important here as it could improve largely the probability 
of a grazing transit, and thus of a larger planet. To inves- 
tigate this point and to assess the validity of the reported 
planetary size, a global analysis of all the data available for 
the HD 80606 system was performed, considering possible 
the presence of trends in the eclipse photometry. This anal- 
ysis presented here aimed also to obtain a reliable estimate 
of the significance of the possible spin-orbit misalignement 
suggested in M09. The confirmation of such a misalignement 
for HD 80606b would be a strong argument for the hypoth- 
esis that the planet owes its peculiar present orbit to the 
Kozai migration mechanism (Wu & Murray 2003). 

The data used in this work and a new photometric re- 
duction of the Spitzer images (L09) are presented in Sec. 2. 
The analysis itself is described in Sec. 3. In Sec. 4, the results 
of the analysis and some of their implications are discussed. 



2 DATA 

The global determination of the system parameters was per- 
formed using the published radial velocities obtained by the 
spectrographs ELODIE (Naef et al. 2001; M09), SOPHIE 
(M09), HRS (Wittenmeyer et al. 2007) and HIRES (Butler 
et al. 2006, L09), the transit photometry obtained in B- 
band by G09 and the one obtained in _R-band by M09 and 
F09. It has to be noticed that the publically available M09 
photometrjQ presents unrealistic error bars, and an extra 
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noise of 0.25 % has to assumed to reproduce the rms of the 
best-fitting solution. 

The occultation observed by L09 brings an important 
constraint on the orbital solution. Instead of using its pub- 
lished timing (L09) as a prior knowledge in the analysis, it 
was chosen to model the L09 Spitzer photometry together 
with the other data. Indeed, an indirect but important con- 
straint on the transit shape comes from the occultation 
shape, the parameters of the two eclipses being of course 
correlated. This is important here, because of the lack of a 
full transit. The Spitzer IRAC occultation photometry pre- 
sented in L09 shows a significant level of correlated noise, 
but this latter is mostly due to the 'per-pixel' photometric 
reduction done in that study to isolate reliably the heat- 
ing of the planet during its periastron passage (D. Deming, 
priv. comm.). A more classical reduction of these data was 
thus performed to better constrain the shape of the occulta- 
tion. Starting from the IRAC images calibrated by the stan- 
dard Spitzer pipeline and delivered to the community as Ba- 
sic Calibrated Data (BCD), the fluxes were converted from 
the Spitzer units of specific intensity (MJy/sr) to photon 
counts, and aperture photometry was obtained for HD 80606 
and HD 80607 in each image using the IRAF/DAOPHOlQ soft- 
ware (Stetson, 1987). The aperture was centered in each 
image by fitting a Gaussian profile on the targets. A mean 
sky background was measured in an annulus extending from 
10 to 20 pixels from the center of the aperture, and sub- 
tracted to the measured fluxes for each image. Each flux 
measurement was compared to the median of the 10 ad- 
jacent images and rejected as outlier if the difference was 
larger than 4 times its theoretical error bar (~ 0.16 %). The 
aperture was adjusted so as to obtain the highest quality on 
the lightcurve of HD 80607, an aperture of 3.5 pixels giving 
the best result. Finally, only 1709 measurements (~ 6.5h) 
encompassing the occultation were kept to make possible 
the modelisation of the combined effect of the heating mod- 
ulation and the IRAC systematic known as the 'ramp' (see 
e.g. Knutson et al. 2008 and references therein) by a simple 
analytical function. Fig. 1 (first panel from the top) shows 
the resulting occultation photometry. 



3 GLOBAL ANALYSIS 

The global analysis was performed with a modified ver- 
sion of the Markov chain Monte Carlo (MCMC) code de- 
scribed in GiUon et al. (2009). MCMC is a Bayesian infer- 
ence method based on stochastic simulations and provides 
the posterior probability distribution of adjusted parame- 
ters for a given model (see e.g. Tegmark et al. 2004, Gre- 
gory 2005, Ford 2006). The Metropolis-Hasting algorithm 
was here introduced into a Gibbs sampler to improve the 
convergence speed of the chains (see e.g. Ford 2006). For 
all the data, a normal distribution was assumed for the er- 
rors. The used model was based on a star and a transit- 
ing planet on a Keplerian orbit about their center of mass. 
More specifically, a classical Keplerian model was used for 

^ IRAF is distributed by the National Optical Astronomy Obser- 
vatory, which is operated by the Association of Universities for 
Research in Astronomy, Inc., under cooperative agreement with 
the National Science Foundation. 



the radial velocities obtained outside the transit in addition 
to a Rossiter-McLaughlin effect model (Gimenez 2006) for 
the SOPHIE measurements obtained during transit, while 
the eclipse model of Mandel & Agol (2002) multiplied by 
a trend/systematic model was used for the different eclipse 
photometric time-series. For both the B-band and _R-band 
time-series, a time-dependent quadratic polynomium was 
used to model the trend. A similar polynomium was also 
enough to model the fiux variation outside the eclipse in 
the IRAC photometry, more sophisticated models did not 
improve the fit. 

The jump parameters of the MCMC analysis were: 
the planet to star radius ratio Rp/R^, the parameter b' = 
acosi/R, {a being the semi-major axis and i the orbital in- 
clination), the transit duration W (first to fourth contact), 
the transit time of minimum light To, the orbital period 
P, the Lagrangian parameters ecosu and esino; (e being 
the eccentricity and uj the argument of periastron), K2 ~ 
K^\~e^P^I'^ [K being the RV orbital semi-amplitude) , 
the products V sin / cos /3 and V sin / sin /3 {V sin I being the 
projected rotational velocity of the star and /3 the spin-orbit 
angle, see Gimenez 2006), 3 trend coefficients for both tran- 
sit time-series and for the IRAC photometry, a different 
systemic radial velocity for each spectrograph, and finally 
the 8 /im occultation depth dF2. For the transit time-series, 
a quadratic limb darkening law was asssumed, with coeffi- 
cients interpolated from Claret's tables (2000; 2004) for both 
filters and for the spectroscopic parameters derived in M09. 
These coefficients were kept fixed in the analysis. 

Due to is Bayesian nature, MCMC allows the easy use 
of priors based on external measurements/constraints. It is 
a major advantage for this analysis, as no complete high- 
quality transit light curve is available to determine unam- 
bigously the transit parameters. The different priors used in 
the analysis were: 

• Vsinl was determined spectroscopically by different 
teams (Naef et al. 2001, Valenti & Fisher 2005). These re- 
sults were combined and the result V sin / = 1.4± 1.0 km s~^ 
was used in a Bayesian penalty added to the merit function 
of the MCMC. 

• Basing on a compilation of spectroscopic data from the 
literature and on their stellar evolution modeling, M09 de- 
rived a value M« = 0.98 ± 0.10 Mq for the stellar mass 
and R, — 0.98 ± 0.07 -Rq for the stellar radius. The result- 
ing stellar density p, = 1.04 ± 0.25 p© was used here in a 
Bayesian penalty added to the merit function. The stellar 
density itself is not a jump parameter in the MCMC, but is 
derived at each step from several jump parameters. On its 
side, the stellar mass was allowed to vary normally within its 
error bar to propagate its uncertainty to the other physical 
parameters. 

• HD 80606 was photometrically monitored from Arizona 
by the MEarth project (Irwin et al. 2009) the night before 
the detection of the transit from Europe, until ~ 2454876.00 
HJD. J. Irwin reported an absence of ingress in the result- 
ing light curve (G. Laughlin's web site oklo.or£|). Using the 
time sampling and scatter presented in the corresponding 
oklo.org post and assuming an absence of correlated noise. 
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these MEarth data were simulated and the resulting time- 
series was added as input data to the analysis, the goal being 
to give a correct statistical weight to this piece of informa- 
tion instead of rejecting a priori any solution with a transit 
beginning before 2454876.00 HJD. 

MCMC is not a global optimization algorithm, in the 
sense that it aims to quantify the uncertainties of a given so- 
lution and not to fully explore the whole model parameters 
hyperspace to locate several well-separated merit function 
minima. To perform in a first stage such a global explo- 
ration, a grid of 1000 different starting points sampling the 
part of the parameters hyperspace close to the published or- 
bital solutions (F09, G09, L09, M09) was set up. For each 
cell of the grid, a short MCMC (1000 iterations) was per- 
formed. From the merit function minima obtained for each 
grid, the existence of a large and continuous minimum area 
was noticed. Outside this area, many chains managed to fit 
well the radial velocities and the partial transits, succeeded 
to predict the occultation at the observed timing, but failed 
to reproduce the duration of the occultation. A few chains 
succeeded a good fit of the whole set of data but leaded 
to an unrealistic stellar density (and thus to a larger merit 
function), illustrating the interest to use as many priors as 
possible to avoid solution degeneracies. 

At this stage, a MCMC of 10^ steps starting from a 
random location within the minimum area was launched. 
The residuals of the deduced best-fitting solution were then 
analyzed to determine the jitter noise of the radial veloc- 
ities and the rod noise of the photometric time-series (see 
Gillon et al. 2009). The deduced values were used to update 
the error bars of the corresponding measurements. 10 new 
chains (of 5 x 10^ iterations each) starting from different lo- 
cations sampling the solution area were then performed. The 
first 20% of each chain were discarded, then the 10 chains 
were combined, their convergence being checked succesfully 
with the Gelman and Rubin statistic test (Gelman & Ru- 
bin 1992). Table 1 gives the resulting best-fitting value, and 
the 1-a (68.3 %) and 3-cr (99.7 %) error bars for each jump 
and physical parameter, while Fig. 1 shows the best-fitting 
model for the eclipses and orbital radial velocity variations. 

As can be seen in Table 1, a grazing transit is rejected 
to a high level of confidence, and the size of the planet is 
thus rather well constrained. Interestingly, a zero value for (3 
is also rejected to a very high level of confidence. It can also 
be noticed that the combined analysis of all the available 
data leads to a significatively refined orbital solution. 



4 DISCUSSION 

The global analysis presented here leads to the secure con- 
clusion that the orbital axis of HD 80606b and the spin axis 
of its host star are misaligned. Such a misalignement was de- 
tected only for the massive planet X0-3b until now (Hebrard 
et al. 2008, Winn et al. 2009). A migration mechanism based 
on planet-planet scattering could explain the misalignment 
of XO-3b (see Winn et al. 2009 and references therein) but 
not the one of the extremely eccentric HD 80606b (Wu & 
Murray 2003). These last authors presented another migra- 
tion scenaxio to explain the peculiar orbit of HD 80606b. 
Under this so-called Kozai migration scenario, HD 80606b 



formed at a few AU from its host star, within a protoplan- 
etary disk very inclined relative to the stellar binary plane. 
HD 80607 would have then induced large eccentricity and 
inclination variations via the Kozai mechanism (see Mazeh 
et al. 1997 and references therein). During the episodes of 
large eccentricity, tidal forces induced by HD 80606 would 
have been large enough to induce an inward migration of 
the planet, finally taking it out of the Kozai cycles. Prom 
then, the orbital evolution of the planet would have been 
totally dominated by tidal circularization. The spin-orbit 
misalignement deduced here is in good agreement with this 
Kozai migration mechanism, as the Kozai cycles should have 
left the planet in a highly inclined orbit. 

Most of the known transiting planets have a very short 
period making possible getting a large amount of good ob- 
servational constraints on their orbital parameters within a 
few lapse of time. For these planets, the observation of a par- 
tial transit can for instance be discarded from the analysis. 
In the case of a 'long' period planet like HD 80606b, every 
piece of data has to be exploited to make possible scientific 
inferences in the first stages of the system characterization. 
The Bayesian method MCMC is very well suited for such 
global analysis as shown by the refinement of the orbital 
parameters obtained for HD 80606b in the present analy- 
sis. Despite the partial nature of the only transit observed 
so far and the possible presence of low-frequency correlated 
noise in the ground-based photometry, the planetary radius 
reveals to be already well constrained by the existing data 
and to be in concord with basic models of irradiated planets 
(e.g. Burrows et al. 2007, Fortney et al. 2007). It is worth 
mentionning that a preliminary MCMC analysis similar to 
the one presented above but performed without the 2 _R- 
band light curves presented in F09 led to a much larger 
uncertainty on the planetary radius {Rp = I.IIq^i Rjup), 
because a good fit of the data could also be obtained with a 
grazing transit solution. It is no more the case when the F09 
data arc added to the analysis, demonstrating the interest of 
a global analysis to maximize the observational constraints. 

Getting a precise full transit light curve of HD 80606b 
and more spectroscopic transit data is of course desirable to 
refine the planetary density and spin-orbit misalignement. 
The measurement of the planetary thermal emission and 
global heating at other wavclcnghts than 8 j^m is also desir- 
able to improve our understanding of the atmospheric prop- 
erties of this fascinating extrasolar planet. 
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Figure 1. From top to bottom: (1) IRAC 8 /im occultation pho- 
tometry (binned per 5 minutes) with the best-fitting model su- 
perimposed. B-band (2) and ij-band (3) transit photometry di- 
vided by the best-fitting trend model and with the best-fitting 
transit model superimposed. The open symbols denote the syn- 
thetic MEarth data. (4) Radial velocities obtained by SOPHIE 
during the transit with 2 models superimposed: the best-fitting 
one (green line) and a model assuming /3 = deg (red dashed 
line). (5) The whole set of radial velocities used in this analysis 
folded using the best-fitting transit ephemeris. 
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Table 1. Best-fitting values and 68.3% and 99.7% confidence 
intervals deduced from the marginalized posterior distribution for 
the jump and physical parameters. 
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